rownames(merge_data) <- merge_data[, 1]
merge_data1 <- merge_data[, -1]  # 删除非数值列
log2_merge_data <- log2(merge_data1 + 1)#对数据进行log2标准化
# 定义分组规则
samples <- colnames(merge_data1)
groups <- ifelse(grepl("^asym", samples), "asym",
                 ifelse(grepl("^ad", samples), "ad",
                        ifelse(grepl("^ctl", samples), "ctl", NA)))


# 初始化存储 P 值的向量，长度为蛋白数量（行数）
p_values <- numeric(nrow(log2_merge_data))
names(p_values) <- rownames(log2_merge_data)

# 遍历每个蛋白（每行）
for (i in 1:nrow(log2_merge_data)) {
  # 提取当前蛋白的数据
  protein_data <- as.numeric(log2_merge_data[i, ])
  
  # 检查每组非零数据点的数量
  group_counts <- table(groups[protein_data > 0])
  valid_groups <- sum(group_counts >= 3)
  # 输出调试信息
  cat("Protein:", rownames(log2_merge_data)[i], "\n")
  cat("Group counts:", group_counts, "\n")
  cat("Valid groups:", valid_groups, "\n")
  
  if (valid_groups < 3) {
    # 若不满足每组至少 3 个非零数据点的条件，P 值设为 NA
    p_values[i] <- NA
  } else {
    # 执行 ANOVA 检验
    aov_result <- aov(protein_data ~ factor(groups))
    # 提取 P 值
    p_values[i] <- summary(aov_result)[[1]]$"Pr(>F)"[1]
  }
}

# 将蛋白名称和对应的 P 值组合成数据框
result_df <- data.frame(Protein = rownames(log2_merge_data), P_Value = p_values)

# 查看结果的前几行
head(result_df)




BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("enrichplot")
BiocManager::install("bitr")
install.packages("bitr")
library(clusterProfiler)
library(ggplot2)
library(org.Hs.eg.db)
library(enrichplot)
library(bitr)

# 在旧版本R中运行：
install.packages("installr")  # 安装工具包
library(installr)
updateR()                     # 自动下载、安装新版本并迁移包

# 加载数据
load("volcano.RData")
# 提取 p < 0.05 的蛋白

library(ggplot2)


load("volcano.RData")
# 提取 p < 0.05 的蛋白
significant_proteins <- prostat[prostat$P < 0.05,]
# 提取蛋白质名称
proteins_names <- significant_proteins$ID
entrez_ids<- bitr(proteins_names, 
                  fromType = "ALIAS",  # 用蛋白质别名
                  toType = "ENTREZID",
                  OrgDb ='org.Hs.eg.db')
# 进行 GO 富集分析
ego <- enrichGO(gene = entrez_ids$ENTREZID,
                OrgDb = 'org.Hs.eg.db',
                ont = "BP",  
                pAdjustMethod = "BH",
                qvalueCutoff = 0.05,
                readable = TRUE)

# 绘制气泡图
dotplot(ego, showCategory = 10) +
  labs(title = "GO Enrichment Analysis (Biological Process)")

